Vector-virus interaction affects viral loads and co-occurrence

Background Vector-borne viral diseases threaten human and wildlife worldwide. Vectors are often viewed as a passive syringe injecting the virus. However, to survive, replicate and spread, viruses must manipulate vector biology. While most vector-borne viral research focuses on vectors transmitting a single virus, in reality, vectors often carry diverse viruses. Yet how viruses affect the vectors remains poorly understood. Here, we focused on the varroa mite (Varroa destructor), an emergent parasite that can carry over 20 honey bee viruses, and has been responsible for colony collapses worldwide, as well as changes in global viral populations. Co-evolution of the varroa and the viral community makes it possible to investigate whether viruses affect vector gene expression and whether these interactions affect viral epidemiology. Results Using a large set of available varroa transcriptomes, we identified how abundances of individual viruses affect the vector’s transcriptional network. We found no evidence of competition between viruses, but rather that some virus abundances are positively correlated. Furthermore, viruses that are found together interact with the vector’s gene co-expression modules in similar ways, suggesting that interactions with the vector affect viral epidemiology. We experimentally validated this observation by silencing candidate genes using RNAi and found that the reduction in varroa gene expression was accompanied by a change in viral load. Conclusions Combined, the meta-transcriptomic analysis and experimental results shed light on the mechanism by which viruses interact with each other and with their vector to shape the disease course. Supplementary Information The online version contains supplementary material available at 10.1186/s12915-022-01463-4.

and in in vivo experiments [18,19], revealing vector immune-related genes whose expression is regulated following viral infection and shows that a crosstalk between vector genes and the virus is imminent for successful viral replication and transmission [20].
Most of the studies have investigated the interaction of a vector with a single virus [20], and a few considered up to two co-occurring viruses [21,22]. However, vectors can usually carry diverse viruses that may interact with each other [23][24][25][26]. Studies on hosts infected by several viruses have shown that coinfection affects viral traits such as virulence and transmission, either by direct virus-virus interaction or indirectly by competition or cooperation [27][28][29]. Therefore, multi-infection may have a profound effect on viral evolution, diversity, and pathogenicity [30,31].
In contrast, the effect of multiple infections on vectors has received far less attention, although is expected to have a considerable effect on the disease epidemiology [25,26]. Therefore, virus-virus interaction is one of the main challenges in virology today [32]. We can expect that as in the case of a host co-infected by several viruses, multi-infection will also have an effect on the vector-virus interaction. Yet, we cannot simply extrapolate from hostvirus studies to vector-viruses, as the two are undergoing different evolutionary processes: while host-pathogen interactions are antagonistic by definition, vector-virus interactions are less definite and can fall anywhere on the continuum between antagonistic and mutualistic [33]. Therefore, the interaction between vector-borne viruses and how they affect their relationship with the vector remains largely unknown.
Varroa destructor (hereafter "varroa") is a parasitic mite which vectors honeybee pathogenic viruses and is routinely co-infected by multiple viruses [34]. The introduction of varroa to the European honey bee (Apis mellifera) has dramatically changed the honeybee viral landscape, leading to worldwide colony collapses threatening global food security [35]. In varroa-free colonies, the viral diversity is high, while viral loads generally remain low [36]. Normally, 1-2 years following varroa introduction to a naive colony, virus diversity and load shifts greatly, giving rise to high levels of a few viral strains as the colony is dwindling to its final collapse [37]. Despite these observations, not much attention has been given to varroavirus interaction, and how it may explain the varroa key role as a vector of honeybee viral disease. In addition, varroa carries over 20 diverse viruses commonly co-occurring ( Fig. 1), including bee-pathogenic viruses, of which some are highly associated with varroa, while others are of unknown pathogenicity to either varroa or bee [38]. Therefore, the varroa-virus system is a unique opportunity to investigate how different viruses interact with each other and with their vector.
Here we studied the interaction of diverse viruses with varroa's transcriptional network, hypothesizing that different viruses have distinct effects. We explored a large set of RNAseq data from the vector aiming (1) to determine if viruses compete or cooperate in their vector and (2) to test whether interactions between individual viruses and the vector correlate with viral epidemiological traits, such as viral load in the vector, and cooccurrence with other viruses. We found no evidence of competition between viruses, but positive correlations in the abundance of specific viruses. Furthermore, we found a strong correlation between viral co-occurrence and the manner in which they interact with the vector's genes, and have experimentally validated some of these predictions. These results show that multi-viral infection can occur not only in the host, but also in the vector, and that these interactions have implications on the virus relationship with its vector and therefore on the disease course.

The viral landscape is heterogeneous across varroa libraries
We found that each of the final 66 varroa RNAseq libraries contains at least three types of viruses (Figs. 1 and 2). The most prevalent viruses belong to the Iflaviridae family, including the bee pathogenic viruses DWVa and DWVb, and the varroa-specific virus, VDV2. On the other hand, three viruses whose presence in varroa was reported before were not detected in any of the libraries (CBPV, VPVL_46, and VPVL_36) (Figs. 1 and 2). Viral load homogeneity across libraries varied greatly between viruses: a few viruses were extremely variable (e.g., DWVa was not detected in a few libraries while reaching < 400,000 TPMs in others), while other viruses (such as ARV-1, ARV-2, and VDV2), were found in somewhat similar loads in most libraries. Interestingly, VDV2 and ARV-1 are the only viruses that were present in all varroa libraries.

No evidence for competition between the different viruses
Multi-infection by several viruses raises the obvious possibility of interactions between them. Interestingly, all significant correlations between viral loads were positive (Pearson correlation followed by Benjamini-Hochberg FDR-correction, P < 0.1) (Fig. 3a).

Varroa modules interact with viral loads
Gene network analysis of the 66 varroa RNAseq libraries clustered the 10,247 varroa genes into 12 modules, each module containing co-expressed genes. The modules were numbered from the largest (module number one which contains 4375 co-expressed genes) to the smallest (module number 12 which contains 37 genes) (Additional file 2). We found significant interactions between specific modules and viruses when correlating module eigengenes to viral loads (transcripts per million, TPM) (Pearson correlation followed by Benjamini-Hochberg FDR-correction, P < 0.1, Fig. 3b). The interaction direction and strength, represented by the correlation coefficient, varied between virusmodule pairs. While VDV2, VDV4, ARV-2, and ARV-1 show positive interaction with some modules (modules 2, 3, 6, and 9), 'bee-viruses' (DWVc and DWVa) show negative interaction with a different module (modules 1 and 9). Interestingly, the only virus that shows both negative and positive interactions with vector's modules is the bee-pathogenic virus, DWVa.

GO terms of interacting modules
The significantly interacting modules consist of genes that are enriched for regulatory-related GO terms such as "Regulation of gene expression" (module 7), "Regulation of metabolic processes" (modules 6 and 9), and immune response-related GO terms, such as "regulation of immunoglobulin secretion" GO:0051023 (module 6). In addition, modules 3 and 6 are enriched for viral and symbiont-related GO terms, such as "viral process" GO:0016032 and "modulation by virus of host process" GO:0019054, and "modulation by symbiont of host cellular process, " GO:0044068. For the list of all significant GO terms in the significantly interacting modules, please see Additional file 3. Of the 17 RNAi-gene homologs recently identified in varroa [39], we found that 13 genes belong to modules interacting with viruses: module 1 (negatively interacting with DWVa and DWVc), module 2 (positively interacting with ARV-1), and module 3 (positively interacting with ARV-2 and ARV-1) (Additional file 4). Correlation between viral loads and varroa modules (eigengenes). c Correlation model between viral load correlations (a), and the distance-matrix of the module-virus correlations (b). In a and b, viruses and modules are ordered according to hierarchical clustering; P values of the Pearson coefficient in a and b are adjusted according to FDR-correction; correlation significance marked by (*) 0.1 < P < 0.01; (**) P < 0.01. For analysis in c, the Mantel test for correlation between two matrices was conducted using 1000 permutations. All viruses' loads are transformed by log10 of the TPM, transcripts per million.

Viruses that are found together interact with the vector's gene co-expression modules in similar ways
Among the module-virus interactions (Fig. 3b), we can detect viruses that share a similar interaction with the same module. For example, VDV2 and VDV4 abundance positively correlated with modules 6 and 9. Interestingly, abundances of these two viruses also significantly correlate with each other (Fig. 3a). Similarly, ARV-2 and ARV-1 also have positively correlated abundances and positive interactions with module 3 ( Fig. 3a and b). Likewise, DWVa and DWVc are positively correlated in abundance and are negatively correlated to modules 1 and 9. We evaluated if this pattern between virus-virus interaction and virus-module interaction can be generalized across all virus-virus-module interactions in our data set. Indeed, we found a significant positive correlation between the two distance matrices (Mantel test for correlation between two distance-matrices [40], Mantel statistic r = 0.44, p < 0.001) (Fig. 3c). In other words, there is a correlation between viral co-occurrence and the manner in which they interact with the vector.

Validating varroa-virus interaction: gene silencing is accompanied by a change in viral load
To experimentally validate the varroa-virus interaction as predicted by the gene-network analysis, we silenced the mite's genes using RNAi and tested its viral load, compared to the control (GFP-treated mites). For silencing, we selected candidate genes based on both the genenetwork analysis and on the literature. From module 7, which interacts with the bee-pathogenic virus DWVa, and contains GO terms of "Regulation of gene expression, " we selected five genes which are highly connected to other genes in the module. In addition, these genes also possessed a relevant annotation, based on literature survey and/or the presence of conserved domain in the predicted coded protein (e.g. immune response-related domains, and genes that were previously reported to interact with host/vector) (for the list of silenced genes, their module-membership and annotation, see Additional file 5). Four of the five tested genes were successfully silenced, i.e., for these genes, a significant reduction in relative gene expression was measured in mites treated with dsRNA, compared to control mites (P < 0.05, Wilcoxon signed-ranks test, followed by FDR-correction, Additional file 6). We then compared the relative viral loads of DWVa, VDV2, and ARV-2 using quantitative PCR (qPCR). DWVa was predicted to negatively interact with the genes in module 7, while VDV2 and ARV-2 were not predicted to interact with the module (Fig. 3b). Of all silenced genes, only mites that were treated with dsRNA of Cuticle protein type 8-like gene (short name: CuP8, accession: LOC111248360) showed a significant change in viral load, compared to control mites. Mites with decreased expression of CuP8 had lower VDV2 and ARV-2 viral loads, compared to control mites (P = 0.02, Wilcoxon signed-ranks test, followed by FDR-correction), while DWVa viral load has also decreased, but not significantly (P = 0.2, Wilcoxon signed-ranks test, followed by FDR-correction) ( Fig. 4a and b).
The dsRNA soaking treatment did not affect mite survival, compared to the control GFP-dsRNA treated mites (Fisher exact test for goodness of fit, p < 0.05) (Additional file 7).

Discussion
Vector-borne viruses rely on another organism for transmission. As the vector plays a key role in viral fitness, we hypothesized that different viruses have distinct effects on a vector, observable as changes in gene expression. This was indeed true for viruses associated Fig. 4 Validating varroa-virus interaction using gene-silencing (RNAi). a Relative expression of CuP8 gene in control (GFP-dsRNA, n=9) and silenced mites (treated with CuP8-dsRNA, n=13). b Relative viral load of VDV2, ARV-2, and DWVa in control and silenced mites. Significant difference between control and silenced mites was tested using the Wilcoxon signed-rank test, followed by FDR-correction. The boxplot shows the interquartile of the data values, the inner thicker line of the box represents the median value, and the dots are potential outliers with varroa mites. Namely, loads of individual viral species were associated with specific changes to the vector's gene expression network (Fig. 3b). Interestingly, co-occurring viruses affected the vector gene expression network in similar ways (Fig. 3c), suggesting a link between viral epidemiological traits (their titer and co-occurrence) and its relationship with the vector. We further propose that vector-virus interactions are evolving rapidly, as we found that closely related viruses have distinct effects on the vector's gene expression. Interestingly, we found no evidence of viral competition within the vector. On the contrary, abundances of some viruses are positively correlated in co-infected vectors, suggesting a potential for viral cooperation. These complex dynamics underscore the role of vector-virus interactions for viral fitness.

Viruses interact with specific modules in the vector's transcriptional network
As we hypothesized, viral presence affects vector gene expression. We found that specific modules in the vector's gene expression network respond to changes in viral titers in a species-specific way (Fig. 3b). This indicates that the host responds to viral presence. Some of these responses may represent antiviral defenses [16,17,19,20] or general stress responses [41] that limit damage to the vector, though it is also possible that viruses trigger additional responses in a way that benefits their spread.
In accordance with this work, we found that the interacting modules involve both specific antiviral and non-specific stress responses. The modules included genes within the RNAi pathway (see Additional file 4), the main arthropod antiviral response [42]. These modules were enriched for infection-specific GO terms such as immune response and viral replication. At the same time, the modules were also enriched for regular cell maintenance GO terms such as cell metabolism, and gene expression regulation. It should be noted that most of the modules responding to viral infection are not obviously related either to stress or antiviral gene expression, and their function vis a vis vector or viral fitness is unclear.
Following the network analysis, we have empirically validated the module-virus interactions using genesilencing experiments of genes with high module connectivity, these are "hub-genes" [43]. While we found that experimentally altering host gene expression does affect viral titer, the interactions were not in the direction the network analysis predicted. Namely, we found that the knockdown of CuP8, a cuticular protein that has been found to assist in plant viral transmission by binding to viruses [44], was accompanied by a significant reduction in viral load of VDV2 and ARV-2, and a non-significant reduction in DWVa (Fig. 4b). Yet, CuP8 was a hub-gene in a module negatively correlated with DWVa (Fig. 3b), suggesting that it's expression should be negatively correlated with that of the virus. The experimental result verifies on one hand the importance of the network hubgenes in the vector-virus interaction, while on the other hand, it illustrates the high complexity of the molecular mechanism underlying the vector response, which may involve other factors such as gene-to-gene transcriptional regulation, and interaction with the environment other microorganisms [25].

A link between viral epidemiological traits and their relationship with the vector
Interestingly, viruses with similar effects on the vector's transcriptional network tended to co-occur (Fig. 3c). This observation can be explained by two non-mutually exclusive explanations. On one hand, a high load of viral RNA and proteins can trigger defense or stress responses by the vector [41,45]. In addition, viruses could trigger these responses as a way to manipulate the vector's ability to spread [46,47]. While the former mechanism is borne out by our data (see discussion of stress and antiviral pathways above), we can make several predictions about the possible manipulation of vectors by the viruses. First, we would expect a more viral species-specific response by the vector, rather than a generalized response caused by increased virus titer. Second, we predict that the effect of the virus on the vector will evolve rapidly, with closely related viruses showing different effects on the vector. We explore these predictions below.

Species-specificity of vector-virus interactions
Viruses indeed show different patterns of interaction with varroa that seem linked to their ecology. For example, some of the more pathogenic bee viruses (DWVa and DWVc) interact with the same module in the host gene expression network (modules 1 and 9, Fig. 3b). These viruses are known to be associated with varroa infestation [34,[48][49][50], and have a positive correlation with each other (Fig. 3a). A different set of modules is affected by viruses that show a high level of expression in varroa, though are detected only at low levels in honey bees (VDV2, ARV-2, and ARV-1 (Fig. 3a)) [51][52][53], suggesting that these viruses may be infecting the mites, rather than using them as a vector. Interestingly, DWVa, a major driver of bee population declines, showed strong interactions with the varroa gene expression network, but modules associated with stress and antiviral responses were largely unaffected, suggesting that this virus may avoid them (Fig. 3b). This could be because the virus does not replicate in varroa [54,55], though it does trigger a few other gene expression responses. These data suggest that viruses interact with varroa in diverse ways and that varroa gene expression is not likely solely driven by their titer.
Our study does not shed light on the mechanisms affecting viral titer in varroa. It could be a direct function of interactions with the varroa gene regulatory network for some viruses, as suggested by our RNAi manipulation (Fig. 4). It is also possible that some viruses may passively accumulate in the mites as a result of feeding on the honey bees. However, mite-feeding behavior may be a trait that viruses could affect making the coevolutionary interactions between bees, mites, and viruses harder to disentangle.

Vector-virus interactions modified by relatively small changes in viral sequence
Virus interaction with its host can drive viral evolution and speciation [56], as host shift is a main force in viral diversification [57,58], and can change the pathogen virulence [59]. Our findings show that this also can occur in vector-virus interactions, as we found that closely related viruses have distinctly different effects on the vector's gene expression. The two most well-known variants of the DWV population, DWVa and DWVb are associated with varroa mites and colony collapses [60]. Varroa infection drives the DWV evolution and the relative prevalence of the two variants [37]. Despite some sequence similarity, DWVb is less associated with the mite's infestation level in the colony, compared to DWVa [61], while on the other hand, DWVb is the only strain with concrete evidence for replicating in the mite [55]. We therefore expected that the two variants, although similar in sequence, will have a different interaction with the mite's gene regulatory networks. Indeed, we found that the two variants show contrasting interactions with the vector modules: while DWVa shows both positive and negative interactions with the mite's modules, DWVb shows no significant interactions with varroa modules (Fig. 3b), which is surprising since DWVb was shown to replicate in the mite [55] and should therefore interact with the vector genes. However, we should note that the analysis can be biased by viruses with high titers, such as DWVa, which extreme high viral loads in many samples may mask other virus-module and virusvirus interactions. Our findings suggest that a relatively small modification in the virus's sequence can lead to a great change in the virus interaction with its vector, and that the vector-virus interactions are continuously and rapidly changing, resulting in a diverse viral community.

Co-infecting viruses do not compete in the vector
The "competitive exclusion principle" states that when competing in the same niche, one species will always suppress the other [62]. This was demonstrated also to happen in vectors, at least in mosquito cell lines co-infected with two viruses [63,64]. In the light of the vector being a limited nourishing source, we could have expected the multi-infecting viruses to compete with each other in the varroa mite. However, we found exactly the opposite, and the only significant intra-viral interactions are positive ones (Fig. 3a). This result, along with the high viral diversity in the varroa mite, supports the model by Leeks et al. [65], which suggests that beneficial multi-viral interactions help to maintain high viral diversity. Still, our findings do not exclude the "short-sighted evolution" model of virulence, which argues that in diverse infections, faster-growing (more virulent) strains are favored because they compete for limited resources [66]. As the virulence of these viruses to varroa mite were not tested so far, we cannot conclude at this point the evolutionary mechanism which led to the observed multi-viral dynamic in the mite.
The positive correlation between some of the viruses may even imply mutualistic interactions, a phenomenon observed before for different strains of West Nile virus co-infecting mosquitoes [67]. A few mechanisms for mutualistic virus-virus relationships were suggested before such as cross-immunity, in multi-viral infection of influenza and other respiratory viral diseases [28,29], and structural protein complementation in measles virus [68,69], and in mutants of Dengue virus infecting mosquito cells [70]. However, the latter is the only example for vector-borne viruses. As viral colonization is the bottleneck for transmission, cooperative interaction between viruses in the vector has direct implications for the viral community dynamics, as they can favor specific viral strains that are not necessarily more virulent in a single, or even double-infected vector. This further emphasizes the need to study multi-viral infections and their molecular mechanism.

Conclusions
Our study contributes to the ongoing investigation of the way viruses interact with their vector, and how this affects the disease course, specifically multi-viral infection, a current major gap in vectorborne viral research [25,26]. Our results imply a link between the virus epidemiological traits and its relationship with the vector. In addition, not only that co-occurring viruses interact with each other, but their abundance may predict the way these viruses will regulate their vector, and potentially, its ability to successfully infect a new host. However, experiments on these viruses' biology and effect on the vector are needed to have the full ecological context of these findings. Hopefully, the gene network pipeline established here can be adopted to other vectorborne diseases, opening the way to study the vector and its associated pathogenic and mutualistic symbionts [71,72].

Methods
In this study, we investigated virus-virus and vector-virus interactions by two approaches: [1] meta-transcriptomic and [2] gene-silencing experiments. In the meta-transcriptomic analysis, we looked at the overall viral landscape in the different libraries, detected vector modules (co-expressing genes using gene-network analysis), and correlated the viruses' abundance to the vector modules. Last, we tested if the virus abundance matrix can predict the virus-vector interaction. In the second step of the study, we experimentally validated the vector-virus interaction on specific virus-gene combinations, selected based on the meta-transcriptomic analysis, by RNAisilencing varroa hub genes and measuring the viral load. All analyses were carried out in the R statistical environment [73]. All meta-transcriptomic analyses are available and reproducible directly from the online supplementary data [74].

Vector-virus interaction, meta-transcriptomic analysis Sequence Read Archive (SRA) data collection
To obtain varroa RNAseq data, we searched for "varroa" term in the SRA databases (NCBI, January 2020) with the following filtering criteria: "RNA" (Source filter), "RNA-seq" (Library Strategy filter) using Illumina technology, and the terms "TRANSCRIPTOMIC" and "VIRAL RNA" (Library Source filter). In total, the filtration yielded 71 libraries of varroa transcriptomes from 11 different studies. The libraries vary significantly in mite species, library preparation method and total library size, number of mites per sample, collected geographical region, sex, physiological stage, and the host species and developmental stage from which the mite was collected (for libraries details please see Additional file 8).

Reads mapping and transcripts quantification
The reads were mapped to both available varroa genome (Vdes_3.0, accession number: GCF_002443255 [75]), and to the genomes of 25 selected viruses (Fig. 1). The alignment and estimation of transcript and virus abundances in transcripts per million (TPM), was carried out using Kallisto [76] (version 0.46.1 with default options). Kallisto pseudo-aligned reads to 35, 659 identified varroa isoforms from 10,247 genes.

Mapping virus reads
The viruses were selected following a survey in the literature and NCBI genome database for viruses related to honey bee and/or varroa mite. Among the hundreds of virus sequences, we included viruses that were previously detected in varroa, in addition to honey bee viruses not detected in varroa before, as a "negative control" (Fig. 1). The final 25 viruses are mostly positive ssRNA viruses, five are negative ssRNA viruses (ARV-1 and ARV-2 (Rhabdoviridae), VOV-1 (Orthomyxoviridae), VDV4 and VDV5 (unclassified)), and another three DNA viruses (one circular dsDNA filamentous, and two ssDNA, found only in varroa mite and not in honey bees [77]). Although many DNA viruses were also found in bees [78], their abundance and importance to varroa/bee health is unknown. Therefore, these sequences were not included in the current study.

Filtering data set
Given the diversity of sources, we wanted to make sure that the input data were as homogeneous as possible. A Principal Component Analysis (PCA) of the different libraries based on varroa gene expression showed that five libraries are obvious outliers (Additional file 9). These were excluded from further analysis: library SRR8867385 from Brattel et al. [79]; libraries SRR5109825 and SRR5109827 from Remnant et al. [80], were deep sequenced for small RNA; and library SRR3927496 from a study by Levin et al. [53], in which a specific virome extraction procedure was implemented prior to library preparation. Last, library SRR533974 by Cornman et al. [81], was made of a pull of 1000 mites, which is exceptionally higher than the number of mites used in most of the studies (a pool of one -five mites). All of these exceptional sample preparation procedures may account for these libraries' deviation from the majority of the data sets in the PCA plot. The remaining 66 libraries are distributed somewhat homogeneously in a subsequent PCA (Additional file 9), and their reads were used for further analyses.

Viral abundance analysis
Viral abundance for all 25 viruses revealed that five viruses were not detected in any of the libraries (CBPV, AFV, ANV, VPVl_46, and VPVL_36) (Fig. 2), and so were not included in further analyses. To check for virus-virus interactions within the mite, we conducted correlation matrix abundance using the Pearson correlation method for all viruses' pairs on the log10 transformed TPM (Fig. 3a).

Weighted gene network co-expression analysis (WGCNA)
We used a network analysis approach to identify groups of genes that share a similar expression pattern across a large set of available varroa transcriptomic data (RNAseq). To construct the gene network, a weighted gene co-expression analysis was carried out using the WGCNA package in R, following the authors' tutorial [43,82]. The WGCNA included 4 main steps: (1) network construction and module detection; (2) correlating modules to external information, the varroa viral load; (3) identifying important genes; and [4] GO term enrichment analysis for varroa modules which interact with the viral load.

(1) Network construction and module detection
Based on the analysis of network topology, for the construction of the network, we set our threshold for merging of modules to 0.25, minimum number of 30 genes per module, and the power β of 12. This power is the lowest for which the scale-free topology fit index curve flattens out upon reaching a high value, in this case, when Rsq reaches 0.886 (Additional file 2). We then performed hierarchical clustering of the genes based on topological overlap (sharing of network neighborhood) to identify groups of genes that co-expressed across libraries, these are the network modules (Additional file 2).

(2) Correlating modules to viral load
To test if the varroa modules interact with the different viruses it carries, we correlated the module eigengenes to the viruses' load (log10TPM). We used the Pearson correlation method and adjusted the p-values for multiple comparisons using the Benjamini-Hochberg method to control the false discovery rate [83] (Fig. 3b).

(3) Identifying hub genes in network modules
Hub genes have high connectivity within a module (Module Membership) [43], and their annotation (based on sequence similarity to homologous genes). The Module Membership is calculated by Pearson correlation of the module eigengene and the gene expression. Genes with high Module Membership and relevant annotation are likely to play a role in the vector-virus interaction and are good candidates for later experimental validation.
(4) GO term enrichment analysis for varroa modules GO terms enrichment analyses for the genes in the significantly interacting modules, were conducted with R package GOstats using the hypergeometric test for the association of categories and genes [84]. The test parameters for each species and each ontology (biological process (BP)) using gene ID from NCBI were as follows: p-value cutoff < 0.05, not conditional, and with detection of over-represented GO terms (testDirection = over).

Correlating virus interaction and virus-varroa interactions
To test if we can predict the virus-varroa interaction given the virus abundance, we used Mantel test for correlation between two distance-matrices [40].

Vector-virus interaction, mites' gene-silencing experiment
In the second step of the study, we experimentally validated specific gene-virus interactions, as predicted by the meta-transcriptomic model. To check if indeed the gene expression is correlated to the varroa viral load, we silenced hub genes with high Module Membership and relevant annotation (see the "Identifying hub genes in network modules" section in the meta-transcriptomic method part). We then checked for the viral load in the silenced mites. For the list of the selected genes for silencing, their accession, and annotation please see Additional file 5.

Mites and honey bee collection
Mites and bees (A. mellifera liguistica) were collected from the same colonies, at the apiary of the Okinawa Institute of Science and Technology (OIST). The hives were not treated against mites and were supplemented with sugar solution and 70% pollen cakes as necessary. Mites were collected from drone and worker pupa of different stages and were kept on bees until soaking, up to five hours from collection.

RNAi silencing of varroa genes
DsRNA preparation For the dsRNA preparation, we first synthesized a T7-promotor attached dsDNA of each of the targeted genes by PCR amplification of cDNA prepared from a pool of 5-10 mites, with specific primers with T7-promotor attached to the 5′ end (see Additional file 10 for the primers information, and section "varroa genes primer design" for details on primer design). For PCR amplification we used Phusion ™ High-Fidelity DNA Polymerase (Thermo Fisher Scientific) in a PTC-200 Peltier Thermal Cycler, MJ Research (BioRad, Toronto, Ontario) with the following steps: an initial denaturation at 98 °C for 30 s followed by 30 cycles of denaturation at 98 °C for 10 s, annealing at 60 °C for 10 s, an extension at 72 °C for 30 s, and a final extension of 72 °C for additional 5 min. We checked that the size of the amplicons matches the expected length, by running 3 μl of the PCR product in 1% agarose gel (135V, 20 min), then verified the sequence by purifying the product and Sanger sequencing on ThermoFisher SeqStudio Genetic Analyzer, using the original reverse primer and BigDye ® Direct Cycle Sequencing kit (Thermo Fisher Scientific), following the manufacturer's instructions. Prior to dsRNA synthesis, we purified the PCR products using a MinElute PCR purification kit (QIAGEN), measured their concentration by Qubit ™ 4 Fluorometer (Life Technologies) with dsDNA HS Assay Kit (Invitrogen), and checked the product size by 4200 Tapestation (Agilent, Tokyo, Japan).
Next, 1200ng of the purified dsDNA with T7-promotor attached was used as a template for the dsRNA synthesis, using MEGAscript ™ T7 Transcription Kit (Thermo Fisher Scientific). We followed the manufacturer's protocol, with slight adjustments. The mix in a total volume of 100μl was incubated overnight at 37 °C, following the addition of TURBO DNase buffer to the reaction, and incubation for another 15 mins at 37 °C. We purified and concentrated the RNA mix using MEGAclear Transcription Clean-Up kit (Thermo Fisher Scientific), and finally measured the RNA concentration by Nanodrop, and checked the product size by 4200 Tapestation (Agilent, Tokyo, Japan). To make sure that the dsRNA effect is specific, we also prepared a negative control dsRNA of a non-target gene, green fluorescent protein (GFP), using pET6Xhn-GFPuv vector (Clontech, Takara) as a template.
Soaking experiment For applying the dsRNA into the mite body, we followed a protocol first developed by Campbell et al. [85] and successfully done by us and others [86][87][88]. We added three mites in a 0.5ml tube containing 20μl dsRNA (2.5μg/μl) in 0.9% NaCl solution. The tubes were kept in 4 °C for 10-15 mins then we checked that all of the mites were soaked, and redipped them into the solution, if needed. The mites were kept in 4 °C overnight (~16 h), dried on a filter paper, then put on a bee pupa (all same age), three mites per bee in a gelatin capsule with perforated top (#1, 0.49ml, HF capsules, Matuya, Japan). Following former studies [88,89], showing optimal silencing effectiveness 48 h post dsRNA treatment, the mites were incubated for 48 h in a controlled, dark environment at 34.5 °C, 60-75% RH, and the pupa was replaced after 24 h. After incubation, each movingviable mite was separated in a 1.5-ml tube, snap-frozen in liquid nitrogen, and kept in −80 °C until RNA extraction. The experiment was replicated in seven experimental batches, between October and November 2020, and each batch included control group mites soaked in GFP-dsRNA of the same concentration and kept under the same conditions as described above. We checked the effect of the dsRNA treatment on mite viability, by comparing the numbers of live and dead dsRNA-treated mites to that of the control mites in each experimental batch (Fisher's exact test for goodness of fit, p < 0.05) (Additional file 7).

RNA extraction and cDNA preparation
Each individual mite was processed following the protocol developed in our lab and described previously [90]. Briefly, each individual mite was crushed in a 1.5 ml tube dipped in liquid nitrogen, then RNA was extracted using a slightly modified TRIzol manufacturer's protocol, with 50% volume of reagents. Total RNA quality and quantity were evaluated using Nanodrop spectrophotometer. Three hundred nanograms of purified RNA was used to synthesize a first-strand cDNA using SuperScript II (Invitrogen) and 1:2 ratio of random hexamer and oligo dT primers following the manufacturer's protocol.

Measuring varroa gene expression and viral load
For both viruses and genes, sets of primers were designed with the NCBI primer design tool (utilizing Primer3 and BLAST), with default parameters and product size set to 100-400bp. Primer's sequence, product length size, and gene IDs or viruses' accession numbers are provided in Additional files 10 and 11 for varroa genes and viruses respectively. The product size of all amplicons was checked by running in 1% agarose gel (135V, 20 min).
Varroa genes primer design For each of the varroaselected genes (Additional file 5), we designed two sets of primers, using the gene mRNA sequence as a template. A first set for dsRNA preparation (as described in the "DsRNA preparation" section), and a second, not overlapping primers-set for gene quantification using qPCR.

Identification of local viruses and RdRp primer design
We targeted the conserved gene of RNAdependent RNA polymerase (RdRp) commonly used for the detection and measurement of different RNA viruses in honey bees and varroa mites [91,92]. As the genome sequence of bee viruses is slightly different across geographical regions [55,81,92], we first wanted to obtain the specific sequence of the strain present in our local mites. For that, for each of the three viruses (VDV2, ARV-2, and DWVa), we amplified and sequenced a wide region of the RdRp gene (amplicon size ~800bp) (see the "DsRNA preparation" section for PCR and sequencing details). To verify the sequence, we nBlasted the product to the nr database (NCBI) [93]. The reverse complementary sequence of viral RdRp was then used as a template to design the qPCR primer sets, as described above. For viruses' amplicon sequences, nBlast results, and primer sequence and position, see Additional file 12.
qPCR To evaluate the relative gene expression and viral load we performed a qPCR using a StepOnePlus Real-Time PCR system (Applied Biosystems Japan, Tokyo, Japan) with TB Green ® Premix Ex Taq ™ II (Tli RNaseH Plus, Takara). The cycling conditions were as follows: 95 °C for 30 s, followed by 40 cycles at 95 °C for 5 s, and 60 °C for 30 s. Data normalization and quantitating was done using StepOnePlus Real-Time system software (Applied Biosystems, Japan), with an automatic threshold set. Small subunit of the ribosomal RNA gene (18S rRNA) of varroa was used as a normalizing gene [94], and a control mite (treated with GFP-dsRNA) from the same experimental batch, used as the normalizing sample. For all qPCR assays a no-template control was included (data not shown). To test for differences in gene expression and viral load of dsRNAtreated mites and control mites, we used a-parametric Wilcoxon signed-ranks test. For all statistical tests, the p-values were adjusted for multiple comparisons using the Benjamini-Hochberg method to control the false discovery rate [83].